clear; clc; close all;
% 读取螺旋桨噪声mp3，计算频谱并找出特征频率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 读取音频文件
[signal, fs] = audioread('example1.wav');  
% 用audioread读入名为*.wav的音频文件
% signal：一个 N×C 的矩阵，N 是采样点数，C 是声道数；fs：采样率（Hz）
if size(signal,2) > 1
    signal = mean(signal, 2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 如果是立体声（C>1），则对左右声道取平均，合成单声道信号
% 预处理：去直流、窗函数
signal = signal - mean(signal); % 去掉直流分量，避免零频（0 Hz）分量过大
N = length(signal);
win = hamming(N);        % 生成长度为 N 的汉明窗，平滑信号两端，减少频谱泄露
xw = signal .* win;     % 将窗函数与信号相乘，得到加窗后的时域信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 计算 FFT
X = fft(xw);               % 对加窗后的信号做快速傅里叶变换，得到复数频谱
P2 = abs(X/N);             % 归一化后取模，得到双边振幅谱
P1 = P2(1:floor(N/2)+1);  % 取单边谱（前半部分），因为实信号的频谱是对称的
P1(2:end-1) = 2*P1(2:end-1);         % 除直流和奈奎斯特频率外，其余分量翻倍，保持能量守恒
f = fs*(0:floor(N/2))/N;   % 构造频率向量f，对应P1中的每个点，范围0–fs/2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 绘图：时域与频域
figure('Position',[100 100 800 600]);% 新建一个 800×600 像素的绘图窗口
% 上半区：绘制去直流后的时域波形
% 横坐标 t 单位秒，纵坐标是幅度
subplot(2,1,1);
t = (0:N-1)/fs;
plot(t, signal);
xlabel('时间 (s)');
ylabel('幅度');
title('原始噪声信号（时域）');
grid on;
% 下半区：绘制幅度谱（dB 量纲），频率范围 0 到 fs/2
subplot(2,1,2);
plot(f, 20*log10(P1));
xlabel('频率 (Hz)');
ylabel('幅度 (dB)');
title('噪声信号频谱');
xlim([0 fs/2]);
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 自动峰值检测：找出几个最大的谱峰
f_low = 50;
f_high = 5000;
idx = find(f >= f_low & f <= f_high);
% 仅在 50–5000 Hz 范围内搜索峰值，排除低频风噪和高频杂音
[peaks, locs] = findpeaks(P1(idx), 'SortStr', 'descend', ...
                         'NPeaks', 5, 'MinPeakDistance', 20);
% 使用 findpeaks ：
%   'SortStr', 'descend' 按峰值高度排序
%   'NPeaks',5 只取前 5 个最高峰
%   'MinPeakDistance',20 保证峰值间至少相隔 20 个频点
peak_freqs = f(idx(locs));
% 将检出的索引 locs 转换为真实频率值
% 在图上标出峰值
hold on;
plot(peak_freqs, 20*log10(peaks), 'rv', 'MarkerFaceColor','r');
for k = 1:numel(peak_freqs)
    text(peak_freqs(k), 20*log10(peaks(k))+3, ...
         sprintf('%.0f Hz', peak_freqs(k)), ...
         'HorizontalAlignment','center');
end
% 在频谱图上用红色倒三角标记这些谱峰，并在其上方标注频率值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输出结果
fprintf('检测到的主要噪声特征频率（Hz）：\n');
disp(peak_freqs');         % 在命令行打印出检测到的几个特征频率列表
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 绘制线性幅度谱
figure(2);
plot(f, P1);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Noise Spectrum (Linear Amplitude)');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
